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Abstract —Big data and the Internet of Things era 
continue to challenge computational systems. Several 
technology solutions such as NoSQL databases have 
been developed to deal with this challenge. In order to 
generate meaningful results from large datasets, analysts 
often nse a graph representation which provides an 
intuitive way to work with the data. Graph vertices 
can represent users and events, and edges can represent 
the relationship between vertices. Graph algorithms 
are used to extract meaningful information from these 
very large graphs. At MIT, the Graphulo initiative 
is an effort to perform graph algorithms directly in 
NoSQL databases such as Apache Accumulo or SclDB, 
which have an Inherently sparse data storage scheme. 
Sparse matrix operations have a history of efficient 
implementations and the Graph Basic Linear Algebra 
Subprogram (GraphBLAS) community has developed a 
set of key kernels that can be nsed to develop efficient 
Unear algebra operations. However, in order to use the 
GraphBLAS kernels, it is important that common graph 
algorithms be recast using the linear algebra building 
blocks. In this article, we look at common classes of 
graph algorithms and recast them into linear algebra 
operations nsing the GraphBLAS building blocks. 

I. Introduction 

The volume, velocity and variety |T] of data being 
collected by today’s systems far outpace the ability to 
provide meaningful results or analytics. A common 
way to represent such large unstructured datasets is 
through a graph representation as they provide an 
intuitive representation of large data sets. In such 
a representation, graph vertices can represent users 
or events and edges can represent the relationship 
between vertices. Many recent efforts have looked 
at the mapping between graphs and linear algebra. 
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In such a mapping, graphs are often represented as 
sparse arrays such as associative arrays or sparse 
matrices using a graph schema. One such effort is 
the Graph Basic Linear Algebra Subprogram (Graph¬ 
BLAS) group which looks to provide a set of kernels 
that can be used to cast graph algorithms as sparse 
linear algebraic operations Q. The abilty to represent 
graph algorithms as linear algebraic operations can be 
greatly beneficial for algorithms scaled for large data 
volume such as those in 0. However, for such 
an initiative to be successful, it is important that the 
proposed linear algebra kernels cover a wide variety of 
graph algorithms that are often used by analysts. This 
article looks at common classes of graph algorithms 
and provides an initial set of graph algorithms recast 
as linear algebraic operations. 

The purpose of our present research effort is to 
enable graph algorithms directly on NoSQL (Not Only 
SQL) databases. Databases such as Apache Accumulo 
or SciDB have become a popular alternative to tra¬ 
ditional relational databases due to their high avail¬ 
ability, partition tolerance and performance. NoSQL 
databases often make use of a key value store or 
store information in triples which are similar to the 
way sparse matrices are stored 0. We see a large 
similarity between our work on performing graph 
algorithms directly on NoSQL databases and research 
on the GraphBLAS specification. The GraphBLAS 
community has proposed an initial set of building 
blocks: 

• SpGEMM; Sparse Generalized Matrix Multiply 

• SpM{Sp}V: Sparse Matrix (Sparse) Vector Mul¬ 
tiply 

• SpEWiseX; Sparse Element-wise Multiplication 

• SpRef; Sparse Reference to a subset 

• SpAsgn; Sparse Assignment to a subset 

• Scale; SpEWiseX with a scalar 

• Apply: Apply a function to each element 


Further, these kernels have been described to work 
on alternate semiring structures such as the tropical 
semiring which replaces traditional algebra with the 
min operator and the traditional multiplication with 
the + operator. This flexibility allows a wide variety of 
graph analytics to be represented using the aforemen¬ 
tioned building blocks. Table summarizes classes of 
graph algorithms that are widely used by the graph 
analytics community. 

With the popularity of NoSQL databases and the 
inherent parallels between the representation of data 
in such databases and sparse arrays, our research 
effort looks at determining how kernels from the 
GraphBLAS specification can be evaluated on NoSQL 
databases. However, in order to ensure that these 
kernels will be able to perform common NoSQL 
database tasks, such as exploration and community 
detection, it is important that the proposed kernels 
are able to express a wide variety of common graph 
analytics. 

A. The Graphulo Initiative 

Graphulo |j^ is an ongoing initiative at the Mas¬ 
sachusetts Institute of Technology that looks at using 
the GraphBLAS kernels on the Apache Accumulo 
database. Accumulo is used for a variety of ap¬ 
plications and has some of the highest published 
performance ||7|. A goal of the Graphulo initiative is 
to use Accumulo server components such as iterators 
to perform graph analytics. In order to provide end 
users with a specification to which they can write their 
algorithms, Graphulo is being written to conform to 
the GraphBLAS specifications. 

B. Paper Outline 

In this article, we present an initial set of common 
graph algorithms recast in the language of sparse 
linear algebra and expressed using the proposed 
GraphBLAS kernels. In Section |I^ we introduce the 
base datatype of NoSQL databases - associative arrays 
- and discuss common schemas used to represent large 
graphs in associative arrays. In Section m we recast 
popular graph algorithms from the Exploration & 
Traversal, Subgraph Detection, Centrality and Com¬ 
munity Detection classes of graph algorithms using 
GraphBLAS kernels. In Section|IV]we discuss the re¬ 
sults, limitations and future work and provide readers 
with an understanding of how these algorithms can 
be implemented on NoSQL databases such as Apache 
Accumulo. We conclude the article in Section Ivl 


11. Associative Arrays and Graph Schemas 

The Graphulo project looks at how graph algo¬ 
rithms can be performed on NoSQL databases. As¬ 
sociative arrays are used as the data type for storing 
and manipulating a large variety of complex datasets. 
In order to represent a dataset using associative arrays, 
we look at a few common schemas that can be used. 

A. Associative Arrays 

Associative arrays are used to describe the rela¬ 
tionship between multidimensional entities using nu¬ 
meric/string keys and numeric/string values. Associa¬ 
tive arrays provide a generalization of sparse matrices. 
Formally, an associative array A is a map from d sets 
of keys Ki x K 2 x ... x Kd to a value set V with a 
semi-ring structure 

A : ATi X ATa X ... X ATd ^ V, 

where (E, ©, 0,0,1) is a semi-ring with 
addition operator ©, multiplication operator ©, 
additive-identity/multiplicative-annihilator 0, and 
multiplicative-identity 1. Furthermore, associative 
arrays have a finite number of non-zero values which 
means their support supp{A) = A“^(F\{0}) is 
finite. 

As a data structure, associative arrays returns a 
value given some number of keys and constitute a 
function between a set of tuples and a value space. 
In practice, every associative array can be created 
from an empty associative array by simply adding and 
subtracting values. With this definition, it is assumed 
that only a finite number of tuples will have values, 
and all other tuples will have a default value of the 
additive-identity/multiplicative-annihilator 0. Further, 
the associative array mapping should support opera¬ 
tions that resemble operations on ordinary vectors and 
matrices such as matrix multiplication. In practice, 
associative arrays support a variety of linear algebraic 
operations such as summation, union, intersection, and 
multiplication. Summation of two associative arrays, 
for example, that do not have any common row or 
column key performs a union of their underlying non¬ 
zero keys. 

Graphulo database tables are exactly described us¬ 
ing the mathematics of associative arrays ij^. In the 
D4M schema, a table in the Accumulo database is an 
associative array. In this context, the primary differ¬ 
ences between associative arrays and sparse matrices 
are: associative array entries always carry their global 
row and column labels while sparse matrices do not. 
Another difference between associative arrays is that 
sparse matrices can have empty rows or columns 


Algorithm Class 

Description 

Algorithm Examples 

Exploration & Traversal 

Algorithms to traverse or search 
vertices 

Depth First Search, Breadth First Search 

Subgraph Detection & Vertex Nomination 

Finding subgraphs or components 
within a graph 

K-Truss subgraph detection. Clique detec¬ 
tion 

Centrality 

Finding important vertices or 
within a graph 

Betweenness Centrality, Eigen Centrality 

Similaiity 

Finding parts of a graph which are 
similar in terms of vertices or edges 

Graph Isomorphism, Jaccard Index, Neigh¬ 
bor Matching 

Community Detection 

Look for communities (areas of 
high connectedness or similaiity) 
within a graph 

Topic Modeling, Non-negative matrix fac¬ 
torization (NMF), Principle Component 
Analysis, Singular Value Decomposition 

Prediction 

Predicting new or missing edges 

Link Prediction, Emerging community de¬ 
tection 

Shortest Path 

Finding the shortest distance be¬ 
tween vertices or sets of vertices 

Floyd Warshall, Bellman Ford, A* Algo¬ 
rithm, Johnson’s Algorithm 


TABLE I: Classes of Graph Algorithms 


while associative arrays do not. For the purposes of 
this algorithmic work associative arrays are encoded 
as sparse matrices. 

B. Graph Schemas 

The promise of big data is the ability to correlate 
diverse and heterogeneous data sources to reduce 
the time to insight. Correlating this data requires 
putting data into a common frame of reference so 
that similar entities can be compared. The associative 
arrays described in the previous subsection can be 
used with a variety of NoSQL databases such as 
Accumulo and require a schema to convert the dense 
arbitrary data into a sparse associative representation. 
Given the variety of data, there are a few commonly 
used graph schemas Q which we discuss below. 

1) Adjacency Matrix: In this schema, data is 
organized as a graph adjacency matrix which can 
represent directed or undirected weighted graphs. 
In this schema, rows and columns of the adjacency 
matrix represents vertices, and values represent 
weighted edges between vertices. Adjacency matrices 
provide a great deal of functionality and are one 
of the more common ways to express graphs 
through matrices. For graph G = (V, E) where 
V = {vi,V2, ...,v„} and E = {ei, 63 ,..., e^}, the 
adjacency matrix A is a n x n matrix where; 

-j = /^ ^ 

y number of self loops, if i = j 

2) Incidence Matrix: The incidence matrix rep¬ 
resentation of a graph can represent multi-hyper- 
weighted as well as directed and multi-partite graphs 
(multiple edges between vertices, multiple vertices per 
edge and multiple partitions). The incidence matrix 
representation is capable of representing complex 


graphs when compared to the adjacency matrix rep¬ 
resentation. In the incidence matrix representation, 
matrix rows correspond to edges, and matrix columns 
represent vertices, with nonzero values in a row indi¬ 
cated vertices associated with the edge. The value at a 
particular row-column pair represents the edge weight 
and sign is often used to represent direction. There are 
many representations for the incidence matrix, and a 
common format is described below. 

For graph G = (V,E) where V = 

{vi, V 2 ,..., v„} and E = {ei, 03 ,..., 0 ^}, the 
incidence matrix E is a m x n matrix where; 

{ + \g\ if G goes into Vj 
— \&i\ if G leaves vj 
0 otherwise 

3) D4M Schema: The D4M 2.0 Schema |§, 
provides a four associative array solution, (Tedge, 
Tedge^, Tdeg, and Traw), that can be used to 
represent complex data. The edge tables, Tedge and 
Tedge^, contain the full semantic information of the 
data set in the rows and columns of the associative 
arrays. From the schema described in Q, a dense 
database can be converted to a sparse representation 
by exploding each data entry into an associative array 
where each unique column-value pair is a column. The 
Tdeg array maintains a count of the degrees of each of 
the columns of T edge, and Traw is used to store the 
raw data. A more thorough description of the schema 
is provided in ||^. Once in sparse matrix form, the 
full machinery of linear algebraic graph processing 
and detection theory can be applied. Linear algebraic 
operations applied on associative arrays organized 
using the D4M schema can have useful results. For 
example, addition of two arrays represents a union, 
and the multiplication of two arrays represents a 
correlation. 












III. Algorithms 

There are many different graph algorithms that can 
be analyzed. In this section, we present an overview 
of our work in representing the classes of graph 
algorithms presented in Table using kernels from 
the GraphBLAS specification. For the work presented 
in this section, we encode associative arrays as sparse 
matrices. 

A. Centrality 

Of the many centrality metrics, there are a few that 
are particularly well-suited to the GraphBLAS frame¬ 
work. Degree centrality, for example, assumes that 
a vertex’s importance is proportional to the number 
of connections it shares. Given an adjacency matrix, 
A, this can easily be computed via a row or column 
reduction, depending on whether in- or out-degree is 
of interest. 

Other centrality metrics are explicitly linear alge¬ 
braic in their formulation. For example, eigenvector 
centrality assumes that each vertex’s centrality is 
proportional to the sum of its neighbors’ centrality 
scores. This is equivalent to scoring each vertex based 
on its corresponding entry in the principal eigenvec¬ 
tor, which can be computed via the power method. 
Starting with Starting with a random positive vector 
Xq with entries between zero and 1, we iteratively 
compute 

^k + l - AX/j; 

until |xf’+iXfe|/(||x/c+i|| 2 ||xfc|| 2 ) is close to 1. 

Related metrics are Katz centrality and PageRank. 
Katz centrality considers the number of /c-hop paths 
to a vertex, for all k, penalizing those with higher 
distances. This is also computed via an iterative 
procedure in which the /cth-order degree vector is 
computed, and added to an accumulator as follows: 

dk+i = Adfe 
Xk-\-l — X]i -f Q. 

where do is a vector of Is and we use the same 
stopping criterion as eigenvector centrality. PageRank 
simulates a random walk on a graph, with the possi¬ 
bility of jumping to an arbitrary vertex. Each vertex 
is then ranked according to the probability of landing 
on it at an arbitrary point in an infinite random walk. 
If the probability of jumping to an arbitrary vertex 
is 0, then this is simply the principal eigenvector of 
A^D~^, where D is a diagonal matrix of vertex out- 
degrees. If the probability of a jump is a, then we 
compute the principal eigenvector of 

yr^'i-NxN + (1 — a)A'''D~^. 


As with eigenvector centrality, this can be done using 
the power method, where multiplication by a matrix 
of Is can be emulated by summing the vector entries 
and creating a new vector where each entry is equal 
to the resulting value. All of these centrality measures 
rely on doing iterative matrix-vector multiplications, 
which fits nicely within the scope of GraphBLAS. 

There has also been work on casting betweenness 
centrality—where a vertex’s importance is based on 
the number of shortest paths that contain it—in linear- 
algebraic operations 0. Other metrics, such as close¬ 
ness centrality, will be the subject of future work. 

B. Subgraph detection and vertex nomination 

Detection of interesting and anomalous subgraphs 
has been a problem of interest for the computer 
science community for many years. Examples of this 
problem space include vertex nomination (ranking 
vertices based on how likely they are to be associated 
with a subset of “cue” vertices) Hg, planted clique 
detection and planted cluster detection m- 

A problem related to planted clique and planted 
cluster detection is computing the truss decomposi¬ 
tion. A A:-truss is a graph in which every edge is part 
of at least k — 2 triangles. Any graph is a 2-truss, and 
any fc-truss in a graph is part of a (fc — l)-truss in 
the same graph. Computing the truss decomposition 
of a graph involves finding the maximal fc-truss for 
all fc > 2. A recent technique for computing the 
truss decomposition GD can be easily converted into 
linear-algebraic operations. Define the support of an 
edge to be the number of triangles of which the edge 
is a member. The algorithm can be summarized as 
follows: 

1) Compute the support for every edge. 

2) If there is no edge with support less than k — 2, 
stop. 

3) Otherwise, remove an edge with support less 
than k — 2, update the supports of its associated 
vertices, and go to|g 

In a more efficient algorithm is proposed that 
considers the edges in order of increasing support. In 
the linear-algebraic form, all edges are considered at 
once, and the appropriate edges are removed simulta¬ 
neously. 

To see the linear-algebraic algorithm, first consider 
the unoriented incidence matrix E. Each row of E 
has a 1 in the column of each associated vertex. To 
get the support for this edge, we need the overlap 
of the neighborhoods of these vertices. If the rows 
of the adjacency matrix A associated with the two 
vertices are summed, this corresponds to the entries 


that are equal to 2. Summing these rows is equivalent 
to multiplying A on the left by the edge’s row in E. 
Therefore, to get the support for each edge, we can 
compute EA, apply to each entry a function that maps 
2 to 1 and all other values to 0, and sum each row of 
the resulting matrix. Note also that 


A = E'^E-diag(E'^E), 


which allows us to recompute EA after edge removal 
without performing the full matrix multiplication. We 
take advantage of this fact in Algorithm Within 
the pseudocode, Xc refers to the complement of x in 
the set of row indices. This algorithm can return the 
full truss decomposition by computing the truss with 
fc = 3 on the full graph, then passing the resulting 
incidence matrix to the algorithm with an incremented 
k. This procedure will continue until the resulting 
incidence matrix is empty. This algorithm can be 
realized using the GraphBLAS kernels SpGEMM, 
SpMV, and Apply. 


Data; The unoriented incidence matrix E, 
integer k 

Result; Incidence matrix of k-truss subgraph E/, 
initialization; 
d = sum(E) 

A = E~^ E — diag{d) 

R = EA 
s={R== 2)1 
X = find{s < k — 2) 
while X is not empty do 
E^ = E{x,:) 

E = E{xc ,;) 

dx = sum(£;a;) 

R = R{xc, ■■) 

R = R — E[EJEx — diag{dx)] 
s={R== 2)1 
X = find{s < k — 2) 

end 

return E 

Algorithm 1: Algorithm to compute fc-truss using 
linear algebra. 1 refers to an array of Is 


As an example of computing the /c-truss using the 
algorithm described, consider the task of finding the 
3-truss of the graph in Fig. [T] 

The incidence matrix for the graph shown in Fig- 



Fig. 1; Example 5-vertex graph 


ure [T] is 

'1 1 0 0 O' 

0 110 0 

1 0 0 1 0 

0 0 110 
10 10 0 

0 10 0 1 


From E, we can compute A using the relation 
A = E"^E — diag{d) to be; 
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To get the support, we first compute 


110 0 0 
0 110 0 
10 0 10 
0 0 110 
10 10 0 
0 10 0 1 
112 11 ' 
2 1111 
112 10 
2 1110 
12 12 0 
1110 1 


0 1110 
10 10 1 
110 10 
10 10 0 
0 10 0 0 


The support is then given by 


s = {R== 2)1 = 


0 0 10 0 
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0 0 10 0 
1 0 0 0 0 
0 10 10 
0 0 0 0 0 
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Since k = 3, x will be the set of edges where the 
support is less than 1, i.e., x = {6} and Xc = 
{1, 2, 3,4, 5}. Thus, R and E will be set to their first 
5 rows, and the update will be computed as follows: 
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The pattern of 2s in R did not change with the removal 
of edge 6, so the support will not change. Therefore, 
the graph represented by the new incidence matrix is 
a 3-tmss. 

C. Similarity 

Computing vertex similarity is important in appli¬ 
cations such as link prediction OD- One common 
method for determining the similarity of two vertices 
is to compute the Jaccard coefficient. This quantity 
measures the overlap of the neighborhoods of two ver¬ 
tices in an unweighted, undirected graph. For vertices 
Vi and Vj where N{v) denotes the neighbors of vertex 
V, the Jaccard coefficient is defined as 

Given the connection vectors (a column or row in the 
adjacency matrix A) for vertices Vi and Vj (denoted 
as Oi and Uj) the numerator and denominator of 
Equation can be expressed as ajaj where we 
replace multiplication with the AND operator in the 
numerator and the OR operator in the denominator. 
This gives us 

Jij = {cLi Aaj)./{ai V Oj) 

T — /A^ 

’Jij — ■f'^AND /-^OR- 

This, however, would involve computing a dense 
matrix, and we are primarily interested in cases where 
this is impossible. Two phenomena can be exploited 
that will help provide an efficient implementation: the 
symmetry of J and sparseness of Since J is 


symmetric, we can compute only the upper triangular 
part and then add the transpose. First we compute the 
upper triangular part of the numerator in the entry 
wise division. The numerator is A^jvjj-,, which in an 
unweighted graph is the same as computing a standard 
matrix multiplication. We can represent A as L -f ( 7 , 
where L is strictly lower triangular and U is strictly 
upper triangular. Since A is symmetric, L = {7^. 
Thus, we have 

= {L + Uf = + LU + UL + U^ 

= iu^y + u'^ + u^u + uu^ 

It can be verified that is strictly upper triangular 
and, therefore is strictly lower triangular. After 

we compute the upper triangular part of A^, we can 
divide each nonzero value by the number of total 
neighbors of the associated vertices. Exploiting these 
properties, we can compute the Jaccard coefficient as 
described in Algorithm]^ The triu operation extracts 
the upper triangular part of the graph, as in MATLAB. 
Algorithm |^can be computed using the GraphBLAS 
kernels SpGEMM, SpMV, and SpEWiseX. Comput¬ 
ing the upper triangular part of a graph can be done 
through a user-defined function that implements the 
Hadamard product. For example, if ® 
triu{A) = A ® 1 where f{i,j) = {A{i,j) : i < 
j, 0 otherwise]. An example of the computation on 
the graph in Fig. [T] is provided in Fig. 

Data: Adjacency matrix A 
Result: Matrix of Jaccard indices J 
initialization; 
d = sum(A) 

U = triu(A) 

A = UU'^ 

Y = U^U 

J = U^+ triu(A) -f triu(y) 

J = J — diag{J) 

for each nonzero entry Jij in J do 

I Jij — Jij / i.di 4“ dj ^ij') 

end 

J = J -p JT 

Algorithm 2: Algorithm to compute Jaccard index 
using linear algebra. 

D. Community Detection 

Community detection is a class of graph algo¬ 
rithms designed to find community structures within 
a graph. Graph communities often contain dense 
internal connections and may possibly overlap with 
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Fig. 2; Computing Jaccard coefficients of the graph in Fig. 0 In line 2, J = + triu (C/C/^) + triu (C/^{7). 

In line 3, J = J./{di + dj — J). Computing J = J + J~^ removes the order dependence. Computation is on 
non-zero entries in each matrix. 


other communities. Real graphs such as social media 
have been shown to exhibit such community struc¬ 
ture on geography, language, age group, etc. Gg. 
The communities may then be used to suggest or 
recommend new information, connections, or even 
products as recommender systems do for popular 
online marketplaces such as Amazon and Google GD- 
One common method used as a basis for such systems 
is topic modeling. Topic modeling is a very popular 
class of algorithms that provides an intuitive look 
into the topics that make up data. As an example, 
consider a set of documents made up of various terms. 
Application of topic modeling can automatically de¬ 
termine a set of topics, the terms that make up a 
topic and the documents that strongly align with these 
topics. Techniques such as topic modeling have gained 
wide usage for automatic summarization, document 
modeling and can provide users with simple and quick 
insight into a dataset. Topic modeling is a general 
field, and a popular technique for topic modeling is 
non-negative matrix factorization GZ). GD- 

Non-negative matrix factorization (NMF) is a class 
of tools used to factorize a given matrix into two 
matrices. Multiplying these two matrices produces 
an approximation of the original matrix. Consider a 


matrix Amxn to be factored into matrices Wmxk and 
Hkxn where m corresponds to the number of rows of 
A, n corresponds to the number of columns in A, and 
k corresponds to the number of topics. Further, NMF 
enforces the constraint that none of these matrices 
contain any negative elements. 

By definition, 

A = W * H. (2) 

In the above factorization, the columns of W can be 
considered a basis for the matrix A with the rows of 
H being the associated weights needed to reconstruct 
A. The property that W and H are nonnegative is 
useful because it can have physical significance (as 
opposed to negative weights or basis elements). One 
way to find the matrices W, H such that A « W * H 
is through an iterative technique such as the algorithm 
presented in Algorithm |3] 

In order to solve the equations in Algorithm 
it is necessary to find a least squares solution to a 
system of linear equations for W and H. One way 
of doing this is by finding the matrix inverse of 
W~^ * W and H * (both are square matrices) and 
multiplying with the right hand side of the equations. 
One method to find the matrix inverse is typically 

















done by techniques such as the Singular Value De¬ 
composition (SVD). However, in order to make use of 
the GraphBLAS kernels, we present an technique used 
by iterative eigenvalue solvers. In such systems, for 
iteration k: X^+i = Xk * {21 — AXk). The algorithm 
used to find the matrix inverse for a square matrix A 
is given in Algorithm]^ 

Data; Matrix A to invert 
Result: X = A-^ 
initialization; 

||^rou>|| mux,; ) 

W^coiW = maxjXX;, 

Xi= A'^/{\\Aro^\\*\\A,ol\\) 
while ||Xt+i - Xt\\F > e do 
I Xt+I = Xt* {2* Inxn — A* Xt) 

end 

Algorithm 4: Matrix inverse through Iteration. At 
each iteration, we check if the value of Xt+i is close 
to the previous iteration estimate of X. 

Using this formulation, computing the inverse of a 
matrix can be done purely using GraphBLAS kernels. 
Combining Algorithms and we can find compute 
the NMF of a matrix A using only GraphBLAS 
kernels. Where {W'^ *W)~^ and {H are de¬ 

termined by using the relation develop in Algorithm]^ 
In fact, computing the NMF of a matrix using Algo¬ 
rithm]^ will require the GraphBLAS SpRef/SpAsgn, 
SpGEMM, Scale, SpEWiseX, and Reduce kernels. 
The outlined algorithm has been tested against a social 
media dataset and provides intuitive results. 

Eor example. Algorithm was applied to a set of 
words collected from the popular social media website 


Data: Incidence Matrix A (size m x n), number 
of topics k 
Result: W and H 
W = random m x k matrix 
while IIA — IF * iT||F > e do 

Solve H ={W'^ * W)-^ *Afor H 
Set elements in iL < 0 to 0 
Solve W'^ = {H* LfT)-! * iL * at for IF 
Set elements in IF < 0 to 0 
end 

Algorithm 5: NME and Inverse through Iteration. 


Twitter. The algorithm was used to determine common 
themes from approximately 20,000 tweets. By setting 
the number of topics to 5, we were able to determine 
words/tweets that fell into 5 different topics. The 
results from this experiment are shown in Eig. 
Erom a graph perspective, this implies that tweets 
corresponding to these tweets from a community. Eor 
topic 1, as an example, this community represents 
users who tweet in the Turkish language. 

IV. Discussion 

The algorithms presented in this paper demonstrate 
several algorithmic capabilities using the initial set of 
GraphBLAS operations, but there are a few inefficien¬ 
cies that could be improved upon with some additional 
functions. In Algorithm]^ for example, when EA is 
computed, it would be more efficient to only consider 
the additions that yield a 2 in the resulting matrix. This 
could be achieved by replacing the + operator with 
a logical AND, but this would violate the semiring 
axioms. Enabling the ability to use linear-algebraic 


Data; Incidence Matrix A (size m x n), number 
of topics k 
Result; Wand H 
initialization; 

W = random m x k matrix 
while IIA — IF * H\\f > e do 

Solve W~^ *W*H =IFT * A for iL 
Set elements in iL < 0 to 0 
Solve H * = H * A'^ for W 

Set elements in IF < 0 to 0 

end 

Algorithm 3: NME through Iteration. At each step 
of the iteration, we check if the Erobenius norm of 
the difference between A and W * H is less than the 
acceptable error. 
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Fig. 3: Application of algorithm to 20k tweets and modeling with 5 topics. Topic 1 represents tweets with 
Turkish words; topic 2 represents tweets related to dating; topic 3 relates to an acoustic guitar competition in 
Atlanta, GA; topic 4 relates to tweets with Spanish words; and topic 5 represents tweets with English words. 
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machinery with data operations that do not conform 
to the rules for semirings may provide substantial 
speedups. 

Algorithm leverages the symmetry of the graph 
to save some of the unnecessary operations, but some 
values under the main diagonal must still be com¬ 
puted in the process. Since it is fairly common to 
work with undirected graphs, providing a version of 
matrix multiplication that exploits the symmetry, only 
stores the upper-triangular part, and only computes 
the upper-triangular part of pairwise statistics, would 
be a welcome contribution to this effort. 

Algorithm computes the NMF of a matrix A 
which can represent the adjacency matrix of a graph. 
However, calculation of the matrix inverse using this 
method can result in dense matrix operations. Since 
the aim of this step is to solve a least squares problem, 
it would be more efficient to implement this using 
a sparse QR factorization or iterative method that 
preserves the sparsity of the problem as much as pos¬ 
sible. We would welcome community involvement in 
building these methods using the GraphBLAS kernels. 

As a next step in the Graphulo effort, we will extend 
the sparse matrix implementations of the algorithms 
discussed in this article to associative arrays. The 
ability to perform the graph algorithms described 
directly on associative arrays will allow us to im¬ 
plement efficient GraphBLAS operations directly on 
Accumulo data structures. In order to make efficient 
implementations, we will use various Accumulo fea¬ 
tures, such as the Accumulo iterator framework, to 
quickly scan Accumulo tables over servers in parallel 
and perform batch operations such as scaling. 

V. Conclusions 

There are a large variety of graph algorithms that 
can be used to solve a diverse set of problems. 
The Graphulo initiative at the Massachusetts Institute 
of Technology is interested in applying the sparse 


linear algebra kernels of the GraphBLAS specification 
to associative arrays which exactly describe NoSQL 
database tables such as those found in the open source 
Apache Accumulo. Current ongoing work includes 
defining efficient implementations of the algorithms 
discussed in this article, extending the classes of sup¬ 
ported algorithms and providing a library that can per¬ 
form basic operations directly in NoSQL databases. 
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